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ABSTRACT 

We present a novel radiation hydrodynamics code, START, which is a smoothed par- 
ticle hydrodynamics (SPH) scheme coupled with accelerated radiative transfer. The 
basic idea for the acceleration of radiative transfer is parallel to the tree algorithm that 
is hitherto used to speed up the gravitational force calculation in an TV-body system. 
It is demonstrated that the radiative transfer calculations can be dramatically accel- 
erated, where the computational time is scaled as N p log N s for N p SPH particles and 
N s radiation sources. Such acceleration allows us to readily include not only numerous 
sources but also scattering photons, even if the total number of radiation sources is 
comparable to that of SPH particles. Here, a test simulation is presented for a multiple 
source problem, where the results with START are compared to those with a radiation 
SPH code without tree-based acceleration. We find that the results agree well with 
each other if we set the tolerance parameter as # cr i t ^1.0, and then it demonstrates 
that START can solve radiative transfer faster without reducing the accuracy. One of 
important applications with START is to solve the transfer of diffuse ionizing photons, 
where each SPH particle is regarded as an emitter. To illustrate the competence of 
START, we simulate the shadowing effect by dense clumps around an ionizing source. 
As a result, it is found that the erosion of shadows by diffuse recombination photons 
can be solved. Such an effect is of great significance to reveal the cosmic reionization 
process. 
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1 INTRODUCTION 

The radiative transfer (RT) in three-dimensional space 
is virtually a six-dimensional problem for a photon dis- 
tribution function in the phase space. So far, various 
RT schemes have been prop osed, some of which are 
coupled with hydrodynamics (|lliev et alj 120061 [2009). In 
a grid-based scheme, the radiative transfer can be re- 
duced to a five- dimensional problem without the signif- 
icant reduction of accuracy [e.g. the ART s cheme that 
is proposed in [ Nakamoto, Umemura, &; Susal ([2001) and 
llliev et al.l ([2006)]. In a smoothed particle hydrodynamics 
(SPH) scheme, if diffuse scattering photons are neglected, 
the computational cost is proportional to N P N S , where N p 
and iV s are the number of SPH particles and that of radi- 
ation sources, respectively. This type of radiation transfer 
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solver can be coupled with the hydrodynamics (e.g. Susa 
2006; llliev et al. 2006, 2009). 

But, in various astrophysical problems, diffuse scatter- 
ing photons play a significant role, for example, in photoion- 
ization of a highly clumpy medium. The treatment of diffuse 
radiation is a hard barrier owing to its high computational 
cost. Susa (2006) proposed Radiation-SPH (RSPH) scheme, 
in which the radiation transfer based on ray-tracing is cou- 
pled with SPH. Since SPH particles are directly used to 
integrate optical depths, high density regions can be auto- 
matically resolved with high accuracy. The code can actually 
treat multiple radiation sources, but the computational cost 
increases in proportion to the number of radiation sources. 
Hence, it is difficult to include numerous radiation sources. 
For the same reason, recombination photons are hard to 
involve, since each SPH particle should be treated as a ra- 
diation source. Thus, in the applicati ons of RSPH, the on- 
the-spot approximation is assumed (jSpitzerl 1 19781 ). where 
recombination photons are supposed to be absorbed on the 
spot, and therefore transfer of diffuse radiation is not solved. 
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The RSPH has been applied to explore the radiative feed- 
back on firs t generation objec t formation (ISusa &; Umemura 
|2004|. |2006l: ISusal 120071. 120081: lHasegawa Umemura fe Susa 
2009: ISusa. Umemura &; Hasegawal 12009). where ionization 
of hydrogen and photo-dissociation of hydrogen molecules 
H2 by ultraviolet (UV) radiation are key physics. But, dif- 
fuse radiation potentially changes the ionization fraction 
near ionization- front (I-front), and therefore alter H2 abun- 
dance. Also, diffuse radiation can erode neutral shadows be- 
hind dense clumps. Such erosion is significant to elucidate 
the cosmic reionization process. To treat the diffuse radia- 
tion properly, the acceleration of radiation transfer solver is 
indispensable. 

In this paper, we present a novel radiation transfer 
solver for an SPH scheme, START, with the tree-based ac- 
celeration of the radiative transfer. The paper is organized as 
follows. In Section 2, we describe the algorithm in detail. In 
Section 3, we test the new scheme by comparing with a pre- 
vious scheme, for the propagation of ionization fronts around 
multiple sources. In Section 4, we demonstrate the impacts 
of diffuse recombination photons by solving the transfer of 
diffuse radiation. Section 5 is devoted to the conclusions. 
Some discussion on the related topics is also given there. 



2 CODE DESCRIPTION 

2.1 Smoothed Particle Hydrodynamics 

The SPH part in our scheme basically follows IMonaghanl 



1992). The density is described as 



Pr 



^mjW(Yij,hi), 



(1) 



where rrij, r^, hi and W are the mass of j-th particle, the 
distance between i-th particle and j-th particle, the smooth- 
ing length of i-th particle, and the kernel function, respec- 
tively. As for the kernel function, we use the standard spline 
form. The equation of motion for each SPH particle i is given 
by 



dvi 
~dt 
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where gi is the gravitational acceleration and P^j) is 
the pressure of i(j)-th particle, n^ is the artificial vis- 
cosity for which we use the Monaghan ty pe viscosity 
(JMonaghan fe Gingoldl Il983l : IMonaghanl 11992 ). W ZJ is the 
symmetrized kernel given by 



Wa 



^(Tij^ + Win^hj) 



(3) 



([Hernquist &; Katzlll989l ). In our scheme, the gravitational 
force on each SP H particle is ca lculated by using Barnes-Hut 
Tree algorithm ([Barnes fe Hutlll986T ). 

As for the equation o f energy, symmetric forms 
have been freque ntly used (|Monaghan &; Gingoldl Il983l ; 
iHernquist fe Katzl ll989T ). However, it is known that sym- 
metric forms sometime s give rise t o negative temperatures 
around a shock front (jBenzJ ll990|) . Hence, we employ an 
asymmetric form for the equation of energy to avoid an 
unphysical behavior around a shock. Such an asymmet- 
ric form has been employed so far by several authours 
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Figure 1. The method for the ray-tracing adopted in RSPH 
scheme is schematically shown. The large red filled circle indicates 
the target particle to which the optical depth from the radiation 
source is evaluated. Smaller filled circles around the target parti- 
cle indicate particles in the neighbor list of the target particle. 



(ISteinmetz fe Mullerl Il993l : lUmemural Il993l : iThacker et all 
[2000), where the energy equation is given by 



dui 
~dt 
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+ H m 3 \^ + 2 Uz - 



Vij-ViWij, (4) 



with the cooling rate Aj and the heating rate Fj . This form 
conserves energy well (Bena ll99Ct IThacker et al.l 12000). In 
our calculations, the error of energy conservation is always 
less than 1 percent. The energy equation is consistently 
solved with the radiative transfer of UV photons and the 
non-equlibrium chemistry. 

2.2 Radiative Transfer 

2.2.1 RSPH ray-tracing 

The ray-tracing algorithm in our scheme is similar to that 
adopted bv ISusal (j_2006), except for the tree-based accelera- 
tion that we newly impleme nt. H e re, we briefly explain the 
RSPH scheme developed bv lSusal (J200d ). 

The steady radiative transfer equation is given by 

dlu 

dr v 

where I u , r u , and S u are the specific intensity, the optical 
depth, and the source function, respectively. The equation 
has the formal solution given by 
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where I u ,q is the specific intensity at r u = 0, and r' v is the 
optical depth at a position along the ray. In order to reduce 
the computationa l cost, so c alled on-the-spot approximation 
is employed ()Spitzerl Il978), where recombination photons 
are assumed to be absorbed on the spot. Therefore, in this 
approximation, the transfer of diffuse radiation is not solved. 
According to the on-the-spot approximation, equation JB} is 
simply reduced to 



Iu(t v ) 
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(7) 



Hence, the specific intensity at a position with r v can be 
evaluated without any iterative method for diffuse radiation. 
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Figure 2. The new ray-tracing method based on oct-tree structure is schematically shown. Each star and circle respectively indicates 
a radiation source and an SPH particle. In the first step, the tree structure for the distributions of radiation sources is constructed. In 
the 2nd step, a list of SPH particles satisfying the same level of conditions (|12[) and (|13|) is produced for each virtual radiation source. 
In the lower right panel, particles in a list of a virtual radiation source are indicated by circles filled with the same color as the owner 
source of the list. In the 3rd step, optical depths from virtual sources to particles in the lists are evaluated in the same manner as Fig. ^ 



In Fig[TJ the method for integrating an optical depth 
in RSPH scheme is schematically shown. The optical depth 
from the radiation source to the target SPH particle can be 
approximately evaluated by 



^target — Tj P ~T~ ^T^ 



(8) 



where r up is the optical depth from the radiation source to 
the upstream SPH particle that is a member of neighbors of 
the target particle and the closest to the light ray (in other 
words, having the smallest angle <j>). At is the optical depth 
between the target particle and the intersection of the light 
ray with a sphere which centers the the radiation source 
and also includes the upstream particle (n4 in Fig[T]). If r up 
is evaluated in advance, we should calculate only At. Here, 
At is evaluated as 



At = aAl 



^target ~T~ Ti\x 



(9) 



where a, A/, ntarget, and n up are the cross section, the dis- 
tance between the intersection and the target particle, the 
number density at the position of the target particle, and 



the number density at the position of the upstream particle, 
respectively. 

In such a method, we have to know the optical depth 
of the upstream particle in advance. If we calculate T up by 
the integration along the light ray, the computational cost is 
roughly proportional to N p ' , where iV p is the total number 
of SPH particles. Hence, the total computational cost is 



J- cal( 



oc N^ /3 N P N B 



(10) 



where iV s is the number of radiation sources. But, if we uti- 
lize the optical depths for SPH particles in order of distance 
from the radiation source, only one SPH particles is used for 
the new optical depth evaluation. Then, the computational 
cost can be alleviated to 

TcaicRSPH OCiV p JV s . (11) 

Nonetheless, this dependence leads to enormous computa- 
tional cost, if we consider numerous radiation sources. Es- 
pecially, to treat the diffuse scattering photons, iV s should 
be equal to iV p , since all SPH particles are emission sources. 
Then, the computational time is proportional to iVp\ Hence, 
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Figure 3. Test 1 — Propagation of ionization fronts around multiple sources in an optically-thick medium, which is assumed to be static 
and uniform. The hydrogen number density is nn = 10 -3 cm -3 and an initial gas temperature is 100 K. The simulation box size is 
132kpc in linear scale. Upper panels show the temperature, where left three panels are the results by RSPH and right three panels are 
the results by START with # cr it — 1- Lower panels show ionization fraction in a slice through the mid-plane of the simulation box. Here, 
the results at three different times t = 0.5t re c, t = 1.0£ rec , and t = 4.0£ rec are shown, where t Tec is the recombination time. 
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Figure 4. Dependence on the tolerance parameter # cr it for Test 
1. In upper panels, the volume fractions of ionized and neutral gas 
are respectively shown by solid and dashed lines, at t = 0.5t rec , 
t = l.Otrec, and t = 4.0trec- Black lines are the results by RSPH, 
while red, green, and blue lines are the results by START with 
#crit = 0.6, # cr it = 1.0, and # cr it = 1.4. Lower panels show the 
volume fraction of temperature. 



to treat diffuse radiation, further acceleration is indispens- 
able. 



until each cell contains only one radiation source. As the 
next step, we reduce the effective number of the radiation 
sources for a target particle i, using the oct-tree structure. If 
a cell of level n is sufficiently distant from the target particle 
z, the radiation sources in the cell are regard as one virtual 
bright source located at the centre of luminosity in the cell. 
Practically, if the following condition is satisfied, the cell is 
regarded to contain one virtual bright source: 



^cell 

d n 



in — 1 

6 cell 



?crit , 



> I 



(12) 



(13) 



where Z™ ell is the size of the n-th level cell and d n is the 
distance from the target particle to the closest edge of the 
cell, and n — 1 indicates the parent cell of the n-th level cell. 
#crit is the tolerance parameter, which regulates the accuracy 
of resulting radiation fields. Of course, such a judgement is 
done not only for n and n—1 levels but also lower level cells, 
such as n— 2, n— 3, ..., as well. Therefore the effective number 
of radiation sources for a target particle can be dramatically 
reduced. In the present scheme, the luminosity of a virtual 
bright source in the n-th level cell is simply given by 
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2.2.2 START: SPH with Tree-based Acceleration of 
Radiative Transfer 

Here, we propose a new ray-tracing method designed to 
solve the radiation transfer for numerous radiation sources. 
The procedure of the new scheme is schematically shown 
in FigEJ For simplicity, the procedure is shown in two- 
dimensional space, although the scheme is implemented in 
three-dimensional space. The new scheme is based on oct- 
tree structure, which is widely used for the gravitational 
force calculations. 

First, we construct the oct-tree structure for the distri- 
butions of radiation sources, where cubic cells are hierarchi- 
cally subdivided into 8 subcells. This procedure continues 



and the position is determined by 

r = — . 



(15) 



where Lj and Yj are the original luminosity and position of 
radiation source j in the cell, respectively. In some situa- 
tions, the evaluation with the SPH kernel as 



L = 



Y, j L J W(r-r J ) 
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Figure 5. Test 2- Propagation of ionization fronts around multiple sources in an optically-thin medium. The same quantities as Fig. [3] 
are shown. 



are other possibilities. Especially, the latter may be a rea- 
sonable choice for the diffuse radiation for which all SPH 
particles are emitters. 

Each virtual radiation source has a list of SPH parti- 
cles which correspond to the same level of tree satisfying 
conditions (|12p and (|13p . The list can be easily produced by 
using the tree structure for the distributions of SPH particles 
that has been already constructed to compute gravitational 
forces. The optical depth towards an SPH particle in the list 
is evaluated in order of distance from the radiation source 
by using equation (|8]). But, if an SPH particle is located 
on the boundary of the domain that the particles in the list 
compose, the particle does not know r up . Therefore, we have 
to integrate the optical depth towards the particle along the 
light ray from the virtual source. 

By this tree-based algorithm, the number of radiation 
sources can be dramatically reduced to the order of log iV s . 
Hence, with the new method, the computational cost is given 
by 



Tcalc, START OC N p logiV s . 



(16) 



Such weak dependence on iV s allows us to solve radiative 
transfer for diffuse photons emitted by all SPH particles as 
well as multiple sources. It is noted that the new ray-tracing 
method accords with RSPH ray-tracing in the case of # C rit = 
or N s = 1. 

The accuracy of the new scheme depends on the value 
of #crit- We can easily guess that the new scheme is quite 
accurate in the optically thick limit, since the contributions 
of distant sources become negligible. In the optically thin 
limit, we can take # C rit similar to that used in the gravity 
calculation, since the radiation flux is proportional to the 
gravity in this limit. For general cases, we discuss what value 
of #crit should be taken in section 3. 



3 COMPARISON BETWEEN RSPH AND 
START 

3.1 Test calculations for the transfer of radiation 
from multiple sources 

As mentioned above, the accuracy and the computational 
time depend on the tolerance parameter # C rit in START. 
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Figure 6. Same as Fig. 3] except for Test 2. 



Hence, we should chose the parameter so as to achieve rea- 
sonable accuracy without vast computational cost. The best 
way to examine the accuracy is to compare simulation re- 
sults with analytic solutions. Unfortunately, when multiple 
radiation sources exist, it is hard to obtain analytic solu- 
tion. Therefore, we compare our results with other results 
simulated with a reliable scheme. Here, we use the results 
by RSPH ray-tracing method, since the RSPH scheme is 
already tested for a standard problem and the ionization 
structure obtained by RSPH is w ell concordant with ana lytic 
solutions as shown in llliev et al.l (2006) and ISusal (|2006h . 

First, we perform tests for numerous sources, under the 
on-the-spot approximation. We calculate the expansion of 
ionized regions around multiple UV sources. In the radia- 
tion hydrodynamic process of primordial gas, not only the 
ionization of hydrogen but also the photo- dissociation of hy- 
drogen molecules H2 is important. Hence, we solve also H2 
chemistry (see Appendix A) and H2 photo-dissociation with 
a shielding function (see Appendix B). We use 128 3 SPH 
particles and 1, 024 radiation sources in all test calculations 
in this section. 
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TIME=1.0t rec 




Figure 7. Test 3 - Expansion of HII regions with hydrodynamics around multiple sources. The initial conditions are the same as those 
in Test 1. Upper panels depict the hydrogen number density, where left three panels are the results by RSPH and right three panels are 
the results by START with # cr it — 1- Middle panels show the temperature. Lower panels show ionization fractions in a slice through the 
mid-plane of the simulation box. Here, the results at three different times t = 0.5t re c, t = 1.0£ r ec, and t = 4.0t re c are shown, where t rec 
is the recombination time. 



3.1.1 Test 1 - Propagation of ionization fronts in an 
optically-thick medium 

As the first test, we consider multiple ionizing sources in an 
neutral, uniform medum with hydrogen number density of 
77-h = 10 -3 cm -3 and an initial gas temperature of 100 K. 
The simulation box size is 132kpc in linear scale. Here, a 
static medium is assumed, and therefore no hydrodynamics 
is solved. The optical depth at the Lyman limit frequency 
initially is ~ 20 for mean particle separation. The radia- 
tion sources are randomly distributed, and simultaneously 
radiate UV radiation when each simulation starts. Each UV 
source has the blackbody spectrum with effective tempera- 
ture of T e ff = 10 5 K and emits ionizing photons per second 
of JV 7 = 5 x 10 48 s _1 . The Stromgren radius is analytically 
given by 



ft 



= (- 



SNy 



(17) 
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where as is the recombination coefficient to all excited levels 
of hydrogen. In the present test, the Stromgren radius for 
each radiation source is 5.4kpc, if the HII regions around 
sources are not overlaped. Simulations are performed until 
the time t — 4.0£ re c, where t rec is the recombination time 
denned by t rec = 1/(tihojb)- 

In Fig[3] the distributions of temperature and ioniza- 
tion fractions are compared between RSPH and START 
with # cr it = 1. In both simulations, individual HII regions 
gradually expand and finally overlap with each other. As 
shown in the figure, the results are almost identical between 
RSPH and START at any evolutionary stage. In order to see 
the dependence on # C rit quantitatively, in Fig. 2] we plot the 
volume fractions of ionized and neutral components (upper 
panels), and the temperature (lower panels) for C vit — 0.6, 
#crit = 1.0, or #crit = 1.4. Here, the results by RSPH are 



also shown. As for ionized and neutral gas fractions, there is 
little disagreement among different # C rit • In the volume frac- 
tion of temperature, there are a slight discrepancy only in 
low temperature regions. The discrepancy increases with in- 
creasing C rit. But, in ionized regions (T > 10 4 K), the results 
by START agree well with RSPH results even for crit = 1.4. 



3.1.2 Test 2 - Propagation of ionization fronts in an 
optically-thin medium 

In the algorithm of START, a test for relatively optically- 
thin medium is important, since the radiation from dis- 
tant sources contribute significantly to ionization structure. 
Therefore, we carry out a test calculation, where the initial 
gas density is rm = 10 -5 cm -3 , which is 100 times lower 
than that in Test 1. The mean optical depth for SPH parti- 
cle separation is ~ 0.2. Each radiation source is assumed to 
have 5 x iV 7 = 10 44 s _1 so that the corresponding Stromgren 
radius is the same as that in Test 1. The positions of sources 
are the same as those in Test 1. 

In FigO the distributions of temperature and ioniza- 
tion fractions are compared between RSPH and START 
with #crit = 1. Here, the early phase (0.5£ rec ), the expanding 
phase (l.Otrec), and the final equilibrium phase (4.0t rec ) are 
shown. Similar to Test 1, each HII region gradually expands 
and finally overlap each other. However, HII regions show 
more bleary structure compared to Fig(3] Such a difference 
originates from the fact that the initial mean optical depth 
is 100 times smaller than that in the Test 1, and each shape 
of ionization front becomes vaguer. 

In Fig. [6] we plot the volume fractions of ionized and 
neutral components (upper panels), and the temperature 
(lower panels) for cr [t = 0.6, C rit = 1.0, or cr [t = 1.4. The 
results by START with # C rit = 0.6 and # C rit = 1.0 are well 
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Figure 8. Dependence on the tolerance parameter cr [ t for Test 
3. In upper panels, the volume fractions of ionized and neutral gas 
are respectively shown by solid and dashed lines, at t = 0.5£ rec , 
t = l.Otrec, and t = 4.0trec- Black lines are the results by RSPH, 
while red, green, and blue lines are the results by START with 
#crit = 0.6, # cr it = 1.0, and # cr it = 1.4. Lower panels show the 
volume fraction of temperature. 



concordant with RSPH results. In the case of # C rit = 1.4, on 
the other hand, a small peak appears at several 10 3 K in the 
volume fraction of the temperature. Hence, # C rit < 1 seems 
to give more reliable results in this case. 



with 



3.1.3 Test 3 - Expansion of HII regions 
hydrodynamics 

Here, we perform a test calculation coupled with hydrody- 
namics. The radiation transfer of UV photons is consistently 
solved with hydrodynamics. The initial gas density and tem- 
perature are set to be the same as those in Test 1. The lumi- 
nosities and positions of radiation sources are also the same 
as those in Test 1. Therefore, the difference between this 
test and Test 1 is solely whether hydrodynamics is coupled 
or not. 

In FigEl the results by START with crit = 1.0 are com- 
pared with the results by RSPH in terms of hydrogen num- 
ber density, temperature, and ionization fractions. As shown 
in this figure, there is no distinct difference between RSPH 
and START, even if we include hydrodynamics. In FigO we 
show the dependence on # C rit in terms of volume fractions of 
ionized and neutral components (upper panels) , and temper- 
ature (lower panels) at three different phases 0.5£ re c, 1.0t rec , 
and 4.0£rec- As shown in the figure, the volume fractions 
for all simulations look well concordant. In Fig [9] we show 
relative differences in temperature (upper panels) and ion- 
ized gas fraction (lower panels) between START and RSPH. 
Here, three START simulations with # C rit = 0.6, # C rit = 1.0, 
and #crit = 1.4 are compared with RSPH at t — 4.0£ rec . The 
distributions are given for the mid-plane of the computa- 
tional box. The relative difference of a physical quantity is 
evaluated as 



AQ = 



IQs 



Qf 



(18) 



where Qstart and Qrsph are the quantities obtained by 
START and by RSPH, respectively. 




Figure 9. Relative differences in the temperature (upper panels) 
and the ionized gas fraction (lower panels) between START and 
RSPH at the epoch t = 4.0£ rec - From left to right column, the 
relative differences between RSPH and START with CY [ t =0.6, 
#crit — 1-0, and # cr it = 1.4 are shown. 




10° 10" 

Nsource 

Figure 10. Average computational time per one step is shown 
as a function of number of radiation sources for SPH particles of 
N p = 128 3 and N p = 64 3 . The computational time of START 
and RSPH are denoted by red and black lines, respectively. As 
for START, the calculations with # crit = 0.6, 0.8, and 1.0 are 
indicated by dashed, solid, and dotted lines, respectively. 



As seen in Fig|9j the relative differences increases ac- 
cording as the tolerance parameter becomes larger. More- 
over, the relative differences in the neutral regions are higher 
than those in the highly ionized region, since the physical 
quantities in such neutral regions are sensitive to the inci- 
dent ionizing flux. Fig[9]shows that if we choose the tolerance 



parameter of C 



1.0, the relative differences become as 



small as < 10 everywhere. 



3.2 Computational time 

In this section, we show how radiative transfer calculations 
are accelerated by the new ray-tracing method. The calcu- 
lations have been done with one Xeon processor (12Gflops). 
The calculation time is measured in Test 3 with varying the 
number of radiation sources. The calculations are done un- 
til t = 1.0t r ec with - 30 - 40 time steps. In FigUHl the 
mean computational time per step is shown as a function 
of the radiation source number. In this figure, the compu- 
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tational time by RSPH is denoted by black line, while that 
by START with crit = 0.6, 0.8, and 1.0 are denoted by 
red dashed, solid, and dotted lines, respectively. The SPH 
particle number is N p = 128 3 or N p = 64 3 . In the case 
with N p = 64 3 , START is always faster than RSPH. This 
actually comes from the reduction of the effective number 
of radiation sources. In the case of N p = 128 3 , START is 
faster than RSPH as long as N s > 10. The calculation with 
N s = 1, 024 by START is roughly 30 times faster than that 
by RSPH, and we find that the results obtained by START 
are well concordant with those by RSPH. The START allows 
the acceleration by N p log N s for a large number of sources, 
as anticipated. However, START is a bit slower than RSPH 
for iV s < 10. In the START scheme, each radiation source 
(cell) has a list of SPH particles, to which optical depths are 
evaluated. When we calculate the optical depths, we also 
have to evaluate optical depths to some SPH particles that 
are not in the list. The time for calculating this additional 
part increases according as the number of SPH particles in- 
creases. In addition, the effective source number is hardly 
reduced in the case of N s ~ 10. Such an increase of com- 
putational time is noticeable for N p = 128 3 run, but not in 
N p = 64 3 run. Anyway, the acceleration by START allows 
us to treat a large number of sources. 



4 RADIATIVE TRANSFER OF DIFFUSE 
RADIATION 

In the previous sections, we have shown the acceleration by 
START for multiple radiation sources. Here, we present the 
effectiveness of START when we include diffuse radiation, 
where all SPH particles should be emitting sources. We con- 
centrate on the photoionization problem again and include 
diffuse recombination photons. 



4.1 Physical process 

To consider the transfer of recombination photons, we adopt 
the recombination coefficient to all bound levels of hydro- 
gen aA (T) , instead of that to all excited levels of hydrogen 
q.-b{T) in the on-the-spot approximation. We solve the trans- 
fer of recombination photons from each SPH particle to all 
other particles by using equation (0) and oct-tree structure. 
The number of recombination photons emitted per unit time 
by an SPH particle j is expressed by 



N ri 



: ai(Tj)n e jn p jVj, 



(19) 



where Tj, n e , %>, and Vj are the temperature, the elec- 
tron number density, the proton number density, and the 
volume of j'-the particle, respectively. ol\ (T) is the recom- 
bination coefficient to the ground state of hydrogen, i.e., 
ai(T) = (xa(T) — qb(T). The emitted energy per one re- 
combination is ~ Jivl + kTj, where h and k are respec- 
tively the Planck constant and the Boltzmann constant. 
The temperature of photoioniz e d gas is typically ~ 10 4 
(e.g., lUmemura fe Ikeuchil Il983 : iThoul fe Weinberg! [l996h , 
and therefore Hul ^> kTj. Thus, the number of recombina- 
tion photons emitted by particle j per unit time and solid 



angle is approximately given by 



Nrecj/^Trdv 





• • • otherwise, 



(20) 
where 8v — kTj/h. As described in Section 2, if a cell is 
far enough from a target particle, SPH particles emitting 
recombination photons in the cell are regarded as one bright 
emitter. The position and the recombination photon number 
of the bright source labeled a are respectively determined 
by 



}_^j A rec, jTj 



AtvSv 



Here 8v is defined by 



59= -r 



k zZj N rec ,jTj 



h J2jN rec ,j 



(21) 



(22) 



(23) 



As for the recombination photon number for virtual sources, 
using another form by the SPH kernel interpolation instead 
of equation (|22[) might be a reasonable way. Using equation 
J7J) instead of equation JB} , the photoionization rate caused 
by a-th radiation source on i-th particle is given by 



kc, = riHi(r, 



•>// 



^L+5z> 



x a v dvd£t. 



(24) 



Here the integration with respect to solid angles is done 
by considering the effect of geometrical dilution. The total 
photoionization rate at r» is given by 



<WO 



^self H - Mother? 



(25) 



where k se u is the contribution from z-th particle itself and 
Mother is the contributions from all emitting particles except 
for i-the particle. k ther is given by 



Mother — / "'Q 



(26) 



Applying equation © and assuming the isotropic radiation 
field around i-th particle, k se u is given by 



kseu = 47m H i / / -^(Jve~ T " +T »dTldv, 

Jur Jo hv 



(27) 



where r v is the optical depth from the edge to the centre of 
the particle, which is given by 



OVTlHI 



3V_ 

4tt 



1/3 



(28) 



In addition, the source function S u can be provided by 



hvot.\ (T)n e n p 




otherwise, 



(29) 



Here the source function can be regarded as almost constant 
near the centre of particle i. As a result, equation (|27[) can 
simply be given by 



h 



elf = / 



VL+Sv 



ai(T) 
8v 



n e n p (l — e Tu ) dv. 



(30) 
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Figure 11. Test 4 ~ Expansion of an HII region in a static 
medium. Upper three panels are the results by the on-the-spot 
approximation, while lower three panels are the results by solv- 
ing the transfer of diffuse recombination photons. The neutral 
fractions in a slice through the mid-plane of the computational 
box are shown at t = 30Myr (left panels), t = lOOMyr (middle 
panels), and t = 500Myr (right panels). 



Figure 12. The temperature distributions in Test 4. 



Note that equation (|25|) accords with ai(T)n e n p , if the 
medium is very opaque everywhere. It means that the 'net' 
recombination coefficient corresponds to as in that case. We 
also evaluate the photoheating rate caused by the recombi- 
nation photons in a similar way. The photoheating due to 
the recombination photons emitted by i-th particle itself, 
and by all other particles are respectively given by 



r S( 

and 









gig) 

8v 



n e n u +h(is - z/l) (l 



dv, (31) 



Tother = ^HlEa/r' " *W* (32) 

xh(v — v\,)e~ Tu > a <j v dvd£l. 

Finally, the total photoheating rate due to the recombina- 
tion radiation for each SPH particle is given by 



4.2 



1 heat — 1 self ~T~ 1 other- 

Test calculations for the transfer of 
recombination photons 



(33) 



In this section, we explore the effect of recombination pho- 
tons by performing some test calculations. All simulations 
in this section are performed with 128 8 SPH particles. The 
transfer of recombination photons is solved by START with 



4.2.1 



:0.i 



Test 4 ~ Expansion of an HII region in a static 
medium 



In this test, we simulate an expansion of an HII region in 
a uniform static medium. The initial physical parameters of 
this test are the same as those of Test 2 in Cosm ologjcal 
Radiative Transfer Comparison Project (|lliev et alj 12006). 
where the hydrogen number density is nu = 10 -3 cm -3 , the 
gas temperature is T — 10 2 K, and the ionization fraction 
is 1.2 x 10 -3 . The computational box is 13.2 kpc in linear 
scale. One radiation source with an effective temperature 
T e ff = 10 5 K and the number of ionizing photons per unit 
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Figure 13. For Test 4, radial profiles of the ionized fraction, the 
neutral fraction (upper row), and the temperature (lower row) 
at t = 30Myr, t = lOOMyr, and t = 500Myr are shown. The 
profiles obtained by three-dimensional radiation transfer START, 
and one-dimensional spherical symmetric radiation transfer are 
indicated by red and black lines, respectively. Solid lines show 
the results by the on-the-spot approximation, while dashed lines 
show those by the full transfer of recombination photons. 



time as N^ 



1 is set up at the centre of the com- 
[0, 0, 0] kpc). The Stromomgren 



putational box ([x, y, z 

radius estimated by equation (fT7|) is 5.4 kpc. The simulation 

is performed until 500 Myr, which roughly corresponds to 

4£ r eo 

In Fig llll we show the neutral fractions in a slice 
through the mid-plane of the computational box at a fast 
expanding phase (t = 30Myr), a slowing-down phase (t = 
lOOMyr), and the final equilibrium phase (t = 500Myr) . 
Fig |12l gives the temperature distributions. In these figures, 
we compare the results by solving the transfer of recombi- 
nation photons with those by on-the-spot approximation. 
Although the difference between these results is not large 
at t — 30Myr ot t — lOOMyr, the difference is distinct at 
later phase at the final equilibrium phase (t = 500Myr). 
The ionized region simulated with the transfer of recombi- 
nation photons is more extended than that simulated with 
the on-the-spot approximation. 

In order to confirm whether the transfer of recombina- 
tion photons is correctly solved, we compare the results ob- 
tained by START with those by a o ne-dimensional spherica l 
symmetric RHD code developed by iKitavama et al.l ((2004) 
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Figure 14. Test 5 - Expansion of an HII region with hydro- 
dynamics. Upper three panels are the results by the on-the-spot 
approximation, while lower three panels are the results by solv- 
ing the transfer of diffuse recombination photons. The neutral 
fractions in a slice through the mid-plane of the computational 
box at t = 30Myr (left panels), t = lOOMyr (middle panels), and 
t = 500Myr (right panels) are shown. 
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Figure 15. The temperature distributions in Test 5. 

in which the transfer of recombination photons is solved. In 
ID RHD simulations, 600 gas shells are used. In Fig |13l we 
show the radial profiles of the ionized fraction, the neutral 
fraction, and the temperature at each evolutionary phase. 
As shown in this figure, the results of the 3D-RHD simu- 
lations and of the ID RHD simulation are well consistent 
with each other. The difference between the profiles in the 
full transfer and the on-the-spot approximation appears in 
moderate regions in which the neutral fractions are x ~ 0.1. 
In these regions, the ionized fractions calculated by solving 
the full transfer are slightly higher than those calculated 
by assuming the on-the-spot approximation. On the other 
hand, near the radiation source, the neutral fractions ob- 
tained by the full transfer are slightly lower, in contrast to 
the moderate ionized regions. This originates from the fact 
that a a is greater than «b- 

4-2.2 Test 5 - Expansion of an HII region with 
hydrodynamics 

This test is similar to Test 4, but here hydrodynamics is con- 
sistently solved with the radiative transfer. The initial phys- 
ical parameter set of the test is the same as those of Test 
5 in Cosmological Radiative Transfer Comparison Project 



Figure 16. The hydrogen number density distributions in Test 
5. 



where the initial gas density and tem- 
j cm~ 3 and 100K. respectively. The size of 
the computational box is 30 kpc in linear scale. A radiation 
source with T eff = 10 5 K and N 1 = 5 x 10 48 s _1 is set at the 
centre of the computational box ([x, y, z] = [0, 0, 0] kpc). 

The neutral fractions, the temperature, and the hydro- 
gen number density in a slice through the mid-plane of the 
computational box are shown in Figs ll4|[T5l and !161 respec- 
tively. In these figures, we show the distributions of these 
quantities at three different evolutionary phases. In the first 
phase, R-type ionization front (I-front) rapidly propagates 
(t = 30Myr). In the second phase, the transition of I-front 
from R-type to D-type occurs (t = lOOMyr), and finally in 
the third phase, D-type I-front preceded by a shock prop- 
agates (t = 500Myr). These figures show that if hydrody- 
namics is coupled, the resultant ionization structure is in a 
quite good agreement with the results by the on-the-spot 
approximation. 

In Fig |17l we compare the results with a 1D-RHD sim- 
ulation. Here, the averaged radial profiles of the ionized 
fraction, the neutral fraction, the temperature, and the hy- 
drogen number density are shown at three different times 
t = lOMyr, t = 200Myr, and t = 500Myr. At every phase, 
the profiles in 3D-RHD simulations are well concordant 
with 1D-RHD results. Similar to the case without hydrody- 
namics, it is found that the recombination photons slightly 
change the profile of ionized fraction near the I-front. On the 
other hand, the profile of the temperature is not strongly 
affected by the recombination photons. As a result, the gas 
dynamics does not change dramatically in this case, even if 
we consider the transfer of recombination photons. 



4.2.3 Test 6 - Shadowing effects by dense clumps 

Here, we explore the shadowing effect behind dense clumps. 
The physical parameters in this test is the same as those in 
Test 4, except that a dense clump with the radius of r c = 
0.56kpc is placed at 0.8 kpc away in the x-direction from the 
box centre ([0.8, 0, 0] kpc). The hydrogen number density of 
the clump is n c = 2.0 x 10 -1 cm -3 , which is 200 times higher 
than that in the surrounding medium of nu = 10 -3 cm -3 . 
The Stromgren radius r s for the density of clump is given 
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Figure 17. For Test 5, radial profiles of the ionized fraction, 
the neutral fraction (top row), the temperature (middle row), 
and the hydrogen number density (bottom row) at t = lOMyr, 
t = 200Myr, and t = 500Myr are shown. The profiles obtained by 
three-dimensional radiation transfer, and one-dimensional spher- 
ical symmetric radiation transfer are indicated by red and black 
lines, respectively. Solid lines show the results by the on-the-spot 
approximation, while dashed lines show those by the full transfer 
of recombination photons. 
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Figure 19. The temperature distributions in Test ! 
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Figure 18. Test 6 - Shadowing effects by a dense clump. The 
hydrogen number density of the clump and surrounding medium 
is respectively n c = 2.0 x 10 _1 cm — 3 and nn = 10 _3 cm — 3 . The 
clump radius is r c = 0.56kpc. Upper three panels are the results 
by the on-the-spot approximation, while lower three panels are 
the results by solving the transfer of diffuse recombination pho- 
tons. The neutral fraction in a slice through the mid-plane of 
the computational box are shown at t = lOOMyr (left panels), 
t = 200Myr (middle panels), and t = 500Myr (right panels). 



by 



2 ' 



(34) 



where F is the incident photon number flux. Assuming T — 
10 4 K and using F at the centre of the clump, the Stromgren 
radius is roughly r s ~ 0.025kpc. Since r c ^> r s , the clump 
is readily shielded from the incident radiation. Thus, it is 



Figure 20. Same as Figjl8l but for the lower density of surround- 
ing medium as nn = 10 — 4 cm -3 . 

expected that the ionization front is trapped in the clump 
and a shadow is created behind the clump. 

In Figs ll8l and [19] we show the resultant ionized frac- 
tion and temperature in a slice through the mid-plane of 
the box. In each figure, the upper panels show the result 
by on-the-spot approximation, while the lower panels show 
that by the transfer of recombination photons. The ioniza- 
tion front is expectedly trapped in the dense clump and a 
shadow is created. Under the on-the-spot approximation, 
only the transfer of UV photons along the radial direction 
from the radiation source is included. On the other hand, 
when we solve the transfer of diffuse recombination pho- 
tons, the recombination photons from the other directions 
photoinize the shadow. We can observe such an erosion of 
shadow in Figs.ll8landll9l But, the erosion is not so drastic 
in this case. It is related to the mean free path of ionizing 
photons. The mean free path of ionizing photons at Lyman 
limit frequency L m f P is given by 



L mfp = 51.4 x 



nH 



pc. 



(35) 



This implies that the initial mean free path in surrounding 
low density regions is approximately one tenth of the size 
of clump. The size of erosion is just corresponding to the 
the mean free path of ionizing photons. That means that in 
lower density of surrounding medium, the erosion is more 
significant. 

In order to confirm the effect of the mean free path, 
we simulate the case in which the number density of sur- 
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Figure 21. The temperature distributions in the case of Fig |20l 
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Figure 22. Ionization of a medium containing small clumps. 
The hydrogen number density and radius of each clump are 
n c = 2.0 x 10 -1 cm — 3 and r c = lOpc, respectively. The num- 
ber density of surrounding medium is nji = 10 -2 cm -3 . Upper 
three panels are the results by the on-the-spot approximation, 
while lower three panels are the results by solving the transfer 
of diffuse recombination photons. The neutral fraction in a slice 
through the mid-plane of the computational box are shown at 
t = 3Myr (left panels), t = lOMyr (middle panels), and t = 20Myr 
(right panels). 



rounding medium is one tenth of that in previous case. In 
this case, the mean free path is roughly comparable with the 
size of the clump, d We show the ionized fraction in a slice 
through the mid-plane of the box in Figj20] The results by 
the on-the-spot approximation and those by the full transfer 
are shown in the upper panels and the lower panels, respec- 
tively. As shown in the figure, the recombination photons 
can gradually ionize the gas behind the clump, in contrast 
to the higher density case. We show the temperature distri- 
butions in Fig J21l Owing to photoheating by the recombina- 
tion photons, the gas behind the clump is heated up. Since 
the number density assumed here roughly corresponds to 
that in intergalactic medium (IGM) at high redshifts, it is 
expected that the recombination photons play an important 
role for the reionization of the Universe. 

We also demonstrate the effect of recombination pho- 
tons for a medium containing multiple small clumps. In this 
simulation, nine dense clumps are set up on the mid-plane 
(z — 0). Each clump has the hydrogen number density of 
n c = 2.0 x 10 -1 cm -3 and the radius of r c = lOpc. The 
number density of surrounding medium is nu = 10 -2 cm -3 . 




Figure 23. Computational time for solving the transfer of diffuse 
recombination photons is shown as a function of the particle (or 
grid cell) number. The computational time with START, RSPH, 
and ART (a grid-based accelerated radiative transfer) are plotted 
by red, black, and blue lines, respectively. As for START, the 
computational time with CY {t = 0.6, 0.8, 1.0, and 1.2 are indicated 
by short-dashed, solid, dotted, and long-dashed lines, respectively. 



These clumps are illuminated by a UV source, which is lo- 
cated at the box centre. The radius of each clump is be 
comparable with the mean free path of ionizing photons. 
Fig [22] shows the ionization structure under the on-the- 
spot approximation (upper three panels), and that result- 
ing from the transfer of recombination photons (lower three 
panels). With the on-the-spot approximation, shadows by 
dense clumps are crearly created. In this case, UV radiation 
never reaches shadowed regions behind the clumps. In con- 
trast, if the transfer of recombination photons is solved, low 
density regions behind all clumps are highly ionized. Con- 
sequently, the ionization structure is dramatically changed, 
compared to the on-the-spot approximation. 



4.3 Computational time for the transfer of diffuse 
photons 

In Fig |231 we show the average computational time per step 
for the simulations in Test 5 as a function of the num- 
ber of SPH particles. The calculations are performed un- 
til t = lOOMyr with ~ 10 - 30 time steps. In this figure, 
RSPH, START, and ART (grid-based accelerated radiative 
transfer) are compared. As mentioned above, in RSPH ray- 
tracing method, the transfer of ionizing photons from each 
SPH particle to all other SPH particles are solved. Therefore, 
the computational time should follow equation (|11[) . On the 
other hand, in START scheme, distant radiation sources are 
regarded as one virtual bright source. Hence, the calculation 
time should follow equation (|16[) . Note that the density dis- 
tribution is nearly uniform at the initial phase, but becomes 
very inhomogeneous at the end of the calculation as shown 
in Fig. 1171 Nonetheless, the measured computational time 
follows equation (|16[) , as seen in Fig. [23] It turns out that 
START is 100-1000 times faster than RSPH in the case of 
N p ~ 10 5 . More dramatic speed-up is expected in simula- 
tions with larger number of particles, although the degree of 
speed-up would be different according to the distributions of 
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particles. As shown in Fig |23l the larger cr [t is, the shorter 
the calculation time is. The computational time is roughly 
scaled to 0~f t . 

Also, the dependence of computational time for START 
is better than that for a grid-based ray-tracing method, 
e.g., a long characteristics method or a short-characteristics 
method. The details are given in Appendix C. ART is a grid- 
based accelerated radiative transfer scheme. The accuracy of 
ART is close to a long-characteristics method, and the com- 
putational cost is similar to a short characteristics method 
(|lliev et alj|200d ). The computational time with ART is ba- 
sically proportional to Nj id . We also show measured cal- 
culation time of ART in Fig. 1231 Actually, the calculation 
time of ART is roughly proportional to Af gr id 5 ' 3 . As shown 
in Fig[23] START is faster than ART, if the total particle 
(grid) number is greater than ~ 10 4 . 



5 CONCLUSIONS AND DISCUSSION 

We have presented a novel accelerated radiation hydrody- 
namics scheme START, which is an SPH scheme with tree- 
based acceleration of radiation transfer. In this scheme, an 
oct-tree structure is utilized to reduce the effective number 
of radiation sources. The computational time by START is 
roughly proportional to N p log N s , where log iV s is the radi- 
ation source number and N p is the SPH particle number. 
We have shown the accuracy of START is almost equivalent 
to that of a previous radiation SPH scheme, if we set the 
tolerance parameter as # C rit ^ 1.0. The method presented 
in this paper provides a powerful tool to solve numerous 
radiation sources. Also, the new scheme allows us to solve 
the transfer of diffuse scattering photons. In this case, the 
computational time is proportional to N p log N p . We have 
simulated the expansion of an HII region around a radiation 
source in an initially uniform medium. Comparing these re- 
sults with ID spherical symmetric RHD, we have confirmed 
that our new method correctly solve the transfer of diffuse 
recombination photons. In addition, we have explored the 
impacts of recombination photons in the case that dense 
clumps are irradiated by a luminous radiation source. As a 
result, we have found that the recombination photons can 
dramatically change the ionization structure especially be- 
hind clumps, if the size of the clump is smaller than the 
mean free path of the recombination photons. 

START is likely to be available for various astrophys- 
ical issues in which numerous radiation sources are re- 
quired. One of them is the simulation of cosmic reioniza- 
tion. Many simul a tions h ave hitherto been w i dely performed 
[e.g.. [Beak et al.1 (I2009T). Illiev et al.1 (I2007T). JMellema et al.1 
(|20od ). ISantos et al.1 J20Q8h . and JThomas et al.1 (|2009h ]. 
With START, it would be possible to trace the reioniza- 
tion process by solving hydrodynamics consistently coupled 
with the transfer of UV photons from numerous sources. 
Furthermore, the diffuse radiation may play an important 
role when we consider the impacts of UV radiation from a 
massive star o n neigh b oring medium (|Ahn &; Shapiro! 
Whalen et al.l 120081; Hasegawa, Umem ura &; Susal 



2007; 



2009; 



Hasegawa, Umemura &; Kitavamal [2009). We will further 
challenge these issues in forth-coming papers. 
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APPENDIX A: NON-EQUILIBRIUM 
CHEMISTRY 

The non-equilibrium chemistry for e - , p, H, H - , H2, and 
HJ~ is implicitly solve d by the chemical network solver in 
iKitavama et al.l (|200l|). Using the evaluated optical depth 
and omitting the suffix i, the photoionization rate of hydro- 
gen for each SPH particle i is given by 



rCic 



£fe< 



(Al) 



where ki on ,a denotes the radiative contribution from radia- 
tion source a on each particle z, which is represented by 



rCi( 



= nm 



r°° r iu, a 

J UJ, J 



exp(- 



hv 



-a u dQdi/, 



(A2) 



where nm, ^l, Iv,a, and a v are the neutral hydrogen number 
density of the particle z, the Lyman limit frequency, the 
intrinsic specific intensity emitted from the radiation source 
a, and the photoionization cross-section at a frequency v. 
Similarly, the photoheating rate for each SPH particle i is 
denoted by 



£r„ 



where 



I VfCX exp(- 



/■oo r 

x (hv — huL)(Tudftdi/. 



(A3) 



(A4) 



Here the integral with respect to solid angles is done by 
taking into account the effect of geometrical dilution. The 
photodissociation rate for each SPH particle i is also given 
by 



k d 



=x> 



(A5) 



where kdis,& is the contribution from radiation source a, 
which is evaluated by using equation (|B1[) . 



APPENDIX B: H 2 PHOTO-DISSOCIATION 

H2 photo-dissociation is regulated by the transfer of Lyman- 
Werner (LW) band lines of H 2 molecules. The transfer o f LW 
band is explored in detail bv lDraine fc Bertoldil (J1996), and 
the self-shielding function is provided. Here, we employ the 
self-shielding function and evaluate the photo-dissociation 
rate as 



/c dis = 1.13 x 10 8 F LW ,o/sh 



N B 



10 14 cm- 2 



where 

fsh(x) 



= U 



3/4 



X < 1 

X > 1 



(Bl) 



(B2) 



Here -Flw,o is the LW flux in absence of the self-shielding 
effect and Nn 2 is the H 2 column density. The H 2 column 
density from a radiation source to the target particle is given 
by 



where A^h 2 ,u P is the column density from the radiation 
source to the upstream particle, and AHu 2 is given by 

'^H 2 ,up + nH 2 ,tar \ 



AiVi 



h 2 = A/ (- 



(B4) 



Here nH 2 , up and nn 2 ,tar are the H 2 number densities at the 
points of the upstream particle and the target particle, re- 
spectively. 



APPENDIX C: 
TRANSFER 



GRID-BASED RADIATIVE 



In the case where all grids emit radiation, the specific in- 
tensity of each grid cell is evaluated by solving equation ((6} 
along characteristics. 

In a long-characteristics method, the specific intensity 
at each grid point are evaluated for all angular directions. 
Moreover, calculations of the order of iV id are needed to 
solve the equation © along each light ray. Therefore, the 
calculation time is roughly proportional to N _ rid x Nq x N^ , 
where N^ , Nq are the number of bins for <f> direction and 6 
direction, respectively. No ~ N<f, should be of the order of 

1 /3 

~ N gr j d , since all the points in the simulation box should be 
irradiated by all grid points. As a result, the computaional 
time is proportional to iVg rid . 

In a short characteristics method, the optical depth is 
calculated by interpolating values of segments lying along 
the light ray. When the transfer of recombination pho- 
tons is involved, the total number of the incident rays 



■2/3 



i nto a computational box is NJja x Nq x N ( 



N, 



■4/3 

grid 



(Nakamo to, Umemura, fc Susal l2QQL). In addition, calcula- 
tion of the order of N _ rid are needed for each incident ray. 

Therefore, the computational time of the short characteris- 

5/3 
tics method is proportional to iV id . 

Authentic Radiative Transfer (ART) scheme is based 
on a long-characteristics method, but only parallel right 
rays are solved. As a result, the accuracy is close to a long- 
characteristics method, and the com putational cost i s simi- 
lar to a short characteristics method (jlliev et al.ll2006h . ART 
scheme allows us to solve the transfer of diffuse radiation. 
Actually, the scheme is applied to evaluate th e escape frac- 
tion of ionizing photons from a high-z galaxy ([Yaiima et al.l 
l2Q09h . In the ART scheme for a simulation of diffuse radi- 
ation, the specific intensity on each segment along a light 
ray is calculated by using equation J6]). The intesity at a 
grid point is obtained by the interpolation from the neigh- 
bouring light rays. The total number of the light rays is 
6 x iVg/ id x N e x N e - iVg/ id . Thus, the calculation time for 

5/3 
this method is proportional to iV id , as similar to that for 

a short characteristics method. 



A^H 2 ,tar = A^H 2 ,up + AiV H2 



(B3) 



